Elastic theory of Normal-Superfiuid Boundary in trapped Fermi Gases 
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By modeling the normal-superfluid boundary in a trapped polarized Fermi gas as an elastic 
membrane, we calculate the atomic density profile. For small trapping anisotropy, we find that the 
superfluid-normal boundary remains approximately elliptical, and has an aspect ratio different from 
that of the trap. For very prolate clouds the boundary becomes distorted into a capsule-like shape. 
We present an analytic explanation of this shape. 



Despite the small numbers of particles (N ~ 10 5 ), the- 
oretical explanations of cold atom experiments almost 
universally involve only bulk properties [l|, Q ■ One ex- 
ception has been experiments on two-component Fermi 
gases in elongated "cigar shaped" traps where in 

addition to bulk physics, one must also include surface ef- 
fects 0, S 0| ■ Here we present a phenomenological model 
of these surface effects, producing a simple explanation 
of the observed density profiles. 

In the experiments of interest, 7 Li atoms are pumped 
into two internal states (f and j) and trapped in an 
anisotropic harmonic trap. Spin relaxation can be ig- 
nored in these experiments, so the number of atoms in 
each state iVf n is conserved. The attractive short-range 
interactions are tuned near unitarity, where the scatter- 
ing length is infinite, yielding a scale-free interaction and 
universal thermodynamics [8| . At low temperature these 
gases form a supcrfluid where "f atoms pair up with J, 
atoms. When the ratio of spins N^/N^ deviates from 
unity, the system phase separates into a central paired 
region surrounded by a predominately polarized region. 
We present a theory of this phase separated gas. 

Our analysis goes beyond the standard local density 
approximation (LDA) which, as described here, can be 
derived from a hydrodynamic theory. Including only 
bulk physics, local hydrodynamic equilibrium requires 
that the pressure in the trap obeys VP = — nW where 
n is the local density and V is the trapping potential. 
Assuming an isothermal trap, this then requires that 
X7p = — W, where /i = + n\)/2 is the average chem- 
ical potential. In the absence of spin-dependent forces, 
the chemical potential difference h = (/ij — /itj,)/2 is inde- 
pendent of space. A consequence of this hydrodynamic 
assumption is that all properties of the gas are unchanged 
as one moves along contours of fixed V. Experiments at 
Rice see a violation of this requirement the domain 
wall separating the superfluid and normal region does 
not follow an elliptical isopotential contour, rather it is 
"squished" axially, forming more of a capsule shape. 

As discussed in 0, [l(| , this distortion is consistent with 
a phenomenological theory where one includes surface 
tension in the domain wall. In 0], this phenomenologi- 
cal theory was explored using a variational ansatz where 
the domain wall was taken to form an ellipse whose as- 
pect ratio was a variational parameter. Haque and Stoof 



[lfl | further explored the shape of the domain wall by 
parameterizing it in cylindrical coordinates (p,z,(f>) as 
a curve obeying (p/R) 1 + (z/Z) 1 — 1 where J,R, and 
Z are determined variationally. Here we optimize the 
shape of the domain wall without restricting its shape in 
any way. We produce an analytic argument which ex- 
plains the shape of the domain wall. The accuracy of 
our approximations are verified via more sophisticated 
numerical calculations. 

Treating the superfluid-normal interface as an elastic 
membrane, hydrodynamic equilibrium requires [Til ] 



2a H = AP, 



(1) 



where oo is the surface tension, H — (l/2)(l/i?i + I/R2) 
is the mean curvature, expressible in terms of the prin- 
cipal radii of curvature R1/2, and AP is the pressure 
difference between the superfluid and normal gas. 

In our calculations, we treat <7o as a spatially uniform 
parameter, which only depends on the average chemical 
potentials of the superfluid and normal gas at the center 
of the trap. In a more sophisticated theory, o-q should 
depend on the local density at the boundary, the pressure 
drop AP, and the temperature. Given that we use an 
approximate equation of state, we feel that our model of 
the bulk properties is not sufficiently accurate to warrant 
including these dependencies whose quantitative forms 
are not all known. 

As introduced by Chevy [12], for the bulk equation of 
state in the superfluid (S) and normal state (N), we use 
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where £s — 



£jv = 1 and £ rs 0.45 is a univer- 



sal parameter. The chemical potentials are [is = A* and 
Pn = A*T = A* + h- This equation of state is exact for 
the zero temperature superfluid where the polarization 
is zero {n-\ = nj.) — it is however only approximate in the 
normal state where it assumes that the local polarization 
is p = (n T - nj.)/(n T + nj.) = 100%. 

Although experiments at MIT [l3| are consistent with 
the local polarization of the zero temperature normal 
state being as low as p — 35%, the Rice experiments Q 
that we are mainly concerned with, find that p > 95%. 
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We are unaware of an explanation for this difference in 
behavior. 

In the superfiuid at unitarity the only microscopic 
energy-scale is the chemical potential, while the trap pro- 
vides a macroscopic length-scale Rtf — ^J'Zpo/rnu! 2 de- 
fined in terms of the central chemical potential po, and 
the trap V (p, z) = (l/2)(muj 2 p 2 + muj 2 z 2 ). We therefore 
introduce a dimcnsionlcss surface tension 



The anisotropy of the trap will be parameterized in terms 
of the dimensionless parameter A = lo z /uj p . A prolate 
cigar-shaped trap has A < 1. 

We parametrize the surface by the function z(p), 
whence the mean curvature has the form H = 
z ppI (l + z p) 2 , where the subscripts refer to derivatives 
with respect to the radial coordinate. 

With this parametrization, Eq. [T] becomes 

2aH = 2a *» = ' [l - (p 2 + AV)] 5/2 (4) 
(1 + zf) C> 

-[l + h-{p 2 +A 2 z 2 )f 2 

where h = h/po- Solving this differential equation, with 
the constraint that the solution is a closed curve, deter- 
mines the shape of the domain wall separating the normal 
and superfiuid regions. 

We first consider the case of an isotropic trap (A = 1), 
where contours of constant pressure drop are circles of 
radius R. In this ansatz, the differential equation reduces 
to an algebraic equation 

2a\ = * (l-R 2 f 2 -(1 + h -R 2 f\ (5) 

which is readily solved graphically. For sufficiently large 
a (at fixed h) there are no solutions to this equation, re- 
sulting in a homogeneous cloud: either all superfiuid or 
all normal. For smaller a there are two solutions: one, 
at smaller radius, representing a local maximum of the 
free energy, and another at larger radius representing the 
local minimum. The local maximum corresponds to the 
critical droplet associated with nucleation of the super- 
fluid phase. The local minimum is shown in Fig.[TJa). As 
one would expect, at fixed po and h, increasing surface 
tension shrinks the radius of the superfiuid region. 

Next, we consider the case of small trapping 
anisotropy, where we find that the full numerical solution 
to Eq. [4] yields a nearly elliptical domain wall. Assum- 
ing an elliptical ansatz, we can calculate the semi-major 
and semi-minor axes a and b. The curvature of an ellipse 
at the axial maximum is H z — b/a 2 , while the curva- 
ture at the radial maximum is H p = (l/2)(a/6 2 + 1/a). 
Substituting these expressions into Eq. dj evaluated at 
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FIG. 1: Depiction of domain walls between superfiuid and 
normal region. For all figures, h equals 0.1. Top Left (a): 
Isotropic trap. Moving inwards from outside, the dimension- 
less surface tension a takes on the values 0, 0.3, 0.6 and 1. The 
radii are normalized with a non-interacting Thomas-Fermi ra- 
dius (Rtf = v/2/xo /muip). Middle left (b): Comparison of 
the domain wall for a = (dashed curve) and a = 0.5 (solid 
curve) for A = u} z /u P = 0.6. Bottom Left (c): Elongated 
trap with A = 0.6 and a — 0.5. The inner dashed line is the 
location of the boundary with the elliptical ansatz. Right (d): 
Cigar shaped trap with A = 0.02 and a = 0.95. Inner dashed 
line: location of boundary from the analytic model. For both 
(c) and (d), the outer solid line is the extent of the atomic 
cloud, and the inner solid line is the shape of the boundary 
obtained from numerical calculations. 



the axial and radial maxima respectively, yields two al- 
gebraic equations that are solved for a and b. The results 
are plotted for h = 0.1, a = 0.5 and A = 0.6. in Fig. rj|b 
and c). The ellipse ansatz indeed works very well even for 
large surface tensions, and moderately anisotropic traps. 
However, further numerical investigation indicates that 
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We explain this shape by writing the free energy of the 
system as 
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FIG. 2: Top (a): Free energy per unit length, S (z) — H(0), 
shown as a function of the radius of the domain wall for var- 
ious values of the axial position z for A = 0.02, a — 0.95 are 
shown. The dashed curve is the free energy for A = and 
a = 0.95. The thick curve plotted for z/Rtf = 6.28, is where 
the free energy minimum occurs at p = p m (z ma x) = 0.62Rtf 
and p = 0. The free energy per unit length H is normalized 

with respect to H = (y^p-) (I?) 5 Po r tf -Bottom (b): A 
quadrant of the numerically determined shape for A = 0.02 
and a = 0.95. The dashed line is the shape predicted from 
the analytic model for the same parameters. The numerical 
maximum occurs at z = 6.56Rtf and the theoretical maxi- 
mum occurs at z = 6.28Rtf- The value of h was set to 0.1 
in all these calculations. 



the numerical solution deviates noticeably from the el- 
lipse ansatz when A < 0.5. 

For large anisotropies A « 1, the boundary becomes 
distinctly non-elliptical, with a pronounced "flattenning" 
on the axial ends. In the extreme case, illustrated in 
Fig- Hid), the domain wall becomes "boxy" with appar- 
ently sharp corners. If one parameterized the boundary 
by writing the radial coordinate as a function p m (z), this 
function appears to be nearly discontinuous. This boxy 
shape has been seen in experimen ts p erformed at Rice 
Q, and in variational calculations [lOj | . 



2tt / dzE(z) 
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(6) 



pdpP N + apr, 




where the dependence of p m and p e on z is implicit. The 
bulk free energy of the system has contributions from 
the superfluid and the normal gas. For every value of 
z, Ps is integrated outward from to the domain wall 
(p m (z)), and P/v is integrated from the domain wall to 
the edge of the trap (p e (z)). The contribution from the 
surface free energy is adA, where the differential area dA 



dp,, 



dz. The angular term 



is expressed as p m W 1 

has been integrated out. 

In the large aspect ratio limit, the slope dp/ dz of the 
boundary should almost everywhere be small. To low- 
est order we neglect this term, and S becomes only a 
function of p m {z) and not its derivative. Thus to mini- 
mize 3>, we simply need to separately minimize S(z) for 
each z. In Fig. |2Ja) we plot 5 as a function of p m (z) for 
various values of z. The near discontinuity observed in 
our numerical calculation is understood by noting that 
for a particular z, S can have multiple minima. Thus 
when the derivative terms are neglected, the p m (z) will 
discontinuously change from a finite value {p m {z max )) to 
zero when one passes the location where these two min- 
ima have the same energy (for example, see the thick 
curve in Fig. Ufa)). This discontinuity is analogous to 
the physics of a first order phase transition. 

We compare the shape predicted by this analytic model 
to the one found in our numerics. As A — > 0, the er- 
ror (Sz^) between the numerically and theoretically ob- 
tained axial maxima (z max ) approaches a value that de- 
pends on a and h, and scales as p m (z ma x)- This scaling 
can be understood by a variational argument. 

Near the axial maximum, one can expand the radial 
coordinate p m (z) as 

Pm (^mai) 



(Sz(°)y 



-(Sz) 



(7) 



In terms of the new variables Sz 



and p 



pm(Zmax) 



the free energy (Eq. O takes the form 



$(,5z (0) ) = 27r<5zW / dSz {E b (Sz) + E s {Sz)) , (8) 
Jo 

where the free energy density B,(Sz) has contributions 
from the bulk fluids (H&) and the domain wall (S s ). The 
bulk energy density is 

5 6 = (pm(z ma x)) 2 / pdpPs - / pdpPN , 

\Jl-(Sz) 2 Jl-(Sz) 2 J 

(9) 
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FIG. 3: (Color Online) Difference between the numerical 
shape of the domain wall and the lowest order analytic 
model (Sz(p)) away from the sharp corners, for various val- 
ues of A at fixed h(= 0.1) and cr{= 0.95). From bottom 
to top, A = 0.02,0.01, and 0.0025. As A -> 0, the differ- 
ence between the numerical and analytically obtained axial 
maximum (Sz m ) = 0.Qlp m {z max ) (= 0.38Rtf). Note that 
p s * z 0> — j ~ 0(1) in agreement with the variational argu- 
ment. The dashed line is the first order perturbative calcu- 
lation of Sz(p) for h = 0.1, A = and a — 0.95, which almost 
entirely accounts for the difference between the analytic and 
numerical results. 



plot 5z(p) for various values of A in Fig. [3] and show that 
after including the first order corrections, the analytic 
calculation is quite close to what we found numerically. 

To summarize, we have determined the shape of the 
superfluid-normal gas domain wall for a range of trap- 
ping anisotropies and surface tensions. For small trap- 
ping anisotropies, the shape is described by a simple el- 
liptical ansatz. For large trapping anisotropies, like those 
in the experiments at Rice 3(, the domain wall is boxy, 
with sharp corners. We have developed a simple model 
to explain this shape, and verified our results with nu- 
merical calculations. 
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merical data. We also acknowledge relevant discussions 
of experimental issues with Randall Hulet, Wenhui Li 
and Wolfgang Ketterle. This work was supported in part 
by the National Science Foundation through grant PHY- 
0456261. 



and the surface energy density is 




where we use the fact that the surface energy density has 

the form p m J 1 + (^^r^J • Upon integrating the free en- 
ergy densities in Eq. [5J we find the free energy has a min- 
imum, whose location depends on the specific choice of a 
and h, and scales as — — r f»OL We have confirmed 
this scaling by comparing the numerical and analytically 
obtained axial maxima for a range of a and h values. 

Furthermore, as A — > 0, the relative error (Sz^ / z max ) 
scales as A. We can investigate the shape of the do- 
main wall in the region near p = by a perturbative 
expansion. For A = 0, we write Sz(p) = z(p) — z max = 
Sz^ + Sz^p 2 + Sz^p 4 , where we have used cylindri- 
cal symmetry to set the coefficients of the odd powers 
to zero. Expanding both sides of Eq. [3] in powers of p, 
one can solve for the coefficients 5z( 4 \ and Sz^ 2 \ We 
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